The Influence of Network Topology on Sound Propagation in Granular Materials 
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Granular materials, whose features range from the particle scale to the force-chain scale to the 
bulk scale, are usually modeled as either particulate or continuum materials. In contrast with 
either of these approaches, network representations are natural for the simultaneous examination 
of microscopic, mesoscopic, and macroscopic features. In this paper, we treat granular materials 
as spatially-embedded networks in which the nodes (particles) are connected by weighted edges 
obtained from contact forces. We test a variety of network measures for their utility in helping 
to describe sound propagation in granular networks and find that network diagnostics can be used 
to probe particle-, curve-, domain-, and system-scale structures in granular media. In particular, 
diagnostics of meso-scale network structure are reproducible across experiments, are correlated with 
sound propagation in this medium, and can be used to identify potentially interesting size scales. 
We also demonstrate that the sensitivity of network diagnostics depends on the phase of sound 
propagation. In the injection phase, the signal propagates systemically, as indicated by correlations 
with the network diagnostic of global efficiency. In the scattering phase, however, the signal is better 
predicted by meso-scale community structure, suggesting that the acoustic signal scatters over local 
geographic neighborhoods. Collectively, our results demonstrate how the force network of a granular 
system is imprinted on transmitted waves. 

PACS numbers: 64.60.aq, 43.25.+y, 81.05.Rm 



During the past 15 years, techniques from areas of 
physics such as statistical mechanics and nonlinear dy- 
namics have been used to make important advances in 
studying networks across myriad disciplines p. Con- 
versely, the perspective of networks can also play impor- 
tant roles in physical problems, as there is a large class 
of heterogeneous systems such as foams, emulsions, and 
granular materials [U [3] for which the connectivity of the 
constituent elements is an important factor in the devi- 
ation of their behavior from continuum models. In fact, 
the discontinuous nature of granular materials led to the 
early idea of a fabric structure governing the anisotropic 
behavior of such materials [4]-[6]. 

We investigate whether studying the rich and complex 
dynamics of granular materials [7 using network analy- 
sis can provide new insights into the underlying physics. 
This treatment is a natural one, because granular mate- 
rials can be represented as spatially-embedded networks 
[8] composed of nodes (particles) and edges (contacts 
between particles) with definite locations in Euclidean 

space [anni. inFig.jij we show a quasi-two-dimensional 

(quasi-2D) granular system composed of photoelastic 
disks that permits the determination of both the con- 
tact network and the interparticle forces. The forces be- 
tween particles in these systems are non- homogeneous, 
and they form a network of chain-like structures that 
span the system (see Fig. [l^). This force chain network 
has the same topology as the contact network but con- 
tains edges that are weighted by the inter-particle forces 



(Fig. [Tp). This is exciting from a networks perspective, 
as it allows us to study the influence of network topology 
on 'network geometry' in a spatially-embedded system. 
From the perspective of granular materials, earlier work 
suggests that force chains provide the main supporting 
structure for static and dynamic loading (TTJ [12] . 

Because signal propagation in granular and heteroge- 
neous materials flT is of considerable importance to nu- 
merous industrial and natural systems, it has been the 
topic of many investigations. A longstanding question 
is how to reconcile the failure of continuum models of 
granular sound propagation pIlHTT] , as such models fail 
to quantitatively describe important heterogeneous and 
nonlinear features of acoustic speed [T8H2I]- The pres- 
ence of force chains has been suggested as a potential 
confounding phenomenon that might underlie the failure 
of previous physical models of sound propagation [18l[22]. 
Ultimately, it would be beneficial to quantify how the 
pressure or strain state of a system is imprinted on trans- 
mitted waves, and to understand how to use these waves 
to accurately detect buried objects or reservoirs of oil. 

An increasing body of work has used tools from ar- 
eas like network science and computational homology to 
obtain insights on the structural properties of granular 
materials [9l [TOl 123 and other continuous media [24] . In- 
deed, a networks perspective provides a valuable comple- 
ment to the standard ways of studying granular materi- 
als. In the present paper, we analyze experimental data 
using network analysis to investigate the role of force- 
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FIG. 1. [Color online] (A) Image of a 2D vertical aggregate of 
photoelastic disks confined in a single layer. The driver posi- 
tion is marked with an arrow. Several particles are embedded 
with a piezoelectric sensor, for which wires are visible. (B) 
The internal stress pattern within the photoelastic particles 
manifests as a network of force chains. (C) Blue lines show a 
weighted graph, which is determined from image processing 
and overlayed on image (B). An edge between two particles 
(nodes) exists if the two particles are in physical contact with 
each other; the forces between particles give the weights of 
the edges. 

weighted contact networks in sound propagation. The 
use of photoelastic particles combined with high-speed 
imaging allows us to gain insight into internal force struc- 
tures and particle-scale sound propagation that are not 
readily available in ordinary granular materials. We find 
that geographic community structure provides a funda- 
mental constraint to sound propagation, illustrating that 
contact topology alone is insufficient to understand signal 
propagation in granular materials. 

EXPERIMENTS 

We perform experiments on a vertical 2D granular sys- 
tem of bidisperse disks confined between two sheets of 
Plexiglass, which have been slightly lubricated with bak- 
ing powder to reduce friction with the container walls. 
The top of the container is open and the particles are con- 
fined exclusively by gravity. The particles are 6.35 mm 



thick and have diameters di = 9 mm and d2 = 11 mm, 
and are cut from Vishay PSM-4 photoelastic material to 
provide measurements of the internal forces. We show 
example images in Fig. [l] These particles have an elastic 
modulus of E = 4 MPa, and they are sufficiently dissipa- 
tive that propagating sound waves experience an approx- 
imately exponential decay as a function of distance from 
the source. The use of such soft, dissipative particles 
differs from previous work [25H29] , where much harder 
particles have been used. Further details about the ap- 
paratus are described in Ref. [22^. 

We excite acoustic waves from the bottom of the sys- 
tem by sending pulses of five 750 Hz sinusoidal waves 
with a voice coil driver attached to a 20 mm wide plat- 
form; maximum particle displacements are on the order 
of 5 /im. To assess the reliability of network diagnos- 
tics, we repeat the experiments for 17 different particle 
configurations, each of which is obtained by manual re- 
arrangement. We restrict our analysis to a region of the 
system that contains just over N = 400 particles. This 
subsystem corresponds to a region in which vertical force 
gradients are minimized due to the Janssen effect [30] . 

We compute particle positions and forces using two 
high-resolution pictures of the static system and one 
high-speed movie that captures the system dynamics. We 
took one static image without the polarizer (see Fig. [TJ^) 
and used it to determine particle positions and contacts 
|31 . We take a second static image using polarizers 
(see Fig. [l^), and we use this together with the contact 
locations to estimate the normal forces at each contact 
using the methodology described in Ref. j22|. We mea- 
sure the amplitude and location of sound propagation 
using a high-speed camera operating at 4000 Hz; the 
camera records 80 frames of data (20 ms) containing both 
the injection of the signal {0 < t < 40) and its dissipation 
(40 < t < 80). For each particle in each frame, we com- 
pute A/(x,7/, t) = I{x^y^t) — I{x,y,to), which measures 
how much the brightness / of the particle changes with 
respect to its unperturbed brightness. In earlier work 
[22] , we used piezoelectric sensors embedded in a subset 
of the particles to determine that AI is proportional to 
the change in stress on that particle. Using AI allows us 
to follow the propagating signal through all particles in 
the measurement region. 

To determine which particles are in contact, we use the 
positions of the particle centers, which are determined 
from the static image of the system using a Hough trans- 
form. If the distance between two particle centers is less 
than 1.05 times the sum of their radii, we treat the par- 
ticles as being in contact. This method overcounts the 
number of true contacts. However, the effect of such over- 
counting is minimized by the fact that false contacts are 
assigned a force value of almost zero when we apply our 
image-processing techniques to find the contact forces. 
Accordingly, they do not contribute to the structure of 
the weighted network. 
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For each experimental run, we construct both an un- 
weighted (binary) and a weighted network, which corre- 
spond respectively to an underlying contact network and 
a force-chain network (see Fig. [T]). In each type of net- 
work, the nodes represent the particles in the system. In 
the binary network A, an edge exists between node i and 
j (i.e., Aij = 1) if node i is in contact with node j; oth- 
erwise, Aij = 0. The weighted network W contains the 
same edges, but each element Wij now has a value that 
is given by an estimate of the normal force fij between 
particles i and j, normalized by the mean force / of all 
contacts: Wij = fij/ f- 



RESULTS 

We assess the global organization of the networks using 
21 candidate diagnostics for A and 8 candidate diagnos- 
tics for W. We define each diagnostic in Appendix A, 
where we also include descriptions to provide intuition 
about what each of them measures, as well as their pos- 
sible physical significance to the granular system that we 
study. We examine the reliability of these diagnostics 
across experimental runs in Appendix B, and we com- 
pare the binary-network diagnostics to those in a null 
model constructed using an ensemble of random geomet- 
ric graphs (RGGs) [32] in Appendix C. We examine 4 di- 
agnostics (clustering coefficient, geodesic node between- 
ness, optimized modularity, and global efficiency) in fur- 
ther detail. Each diagnostic can be defined for both bi- 
nary and weighted networks, and each is helpful for ob- 
taining insights into a particular type of spatial structure 
in the system: particles (cluster coefficient), curves (be- 
tweenness), meso-scale domains (via community struc- 
ture determined from modularity optimization), and the 
entire system (global efficiency). We describe our results 
in the sections below. 
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FIG. 2. [Color online] Example distributions of several net- 
work diagnostics for a sample granular packing. The net- 
work characteristics that we examine include (A) global effi- 
ciency Eu,^ (B) community structure, which we visualize us- 
ing the quantity X'^{z + 5), where X is the community label, 
(C) geodesic node betweenness B^^ and (D) clustering coeffi- 
cient Cw This figure illustrates their respective sensitivities 
to system-scale, domain-scale, curve-scale, and particle-scale 
structure. The quantity X'^{z + 5) allows us to visualize both 
the community label {X) and the intra- community strength 
z-score (2;, see Eq.[3| simultaneously; we chose the constant 5 
purely for visual clarity. 



it is sensitive to curve-like structures in the network. Fi- 
nally, we find that clustering coefficient [see Eqs. 13 



and 24)] is sensitive to particle-scale features of the net- 
work. In a later section, we report how each of these 
network diagnostics correlates with sound propagation 
through the granular material. 



A key advantage of using network tools to study granu- 
lar materials is that different network diagnostics (which 
we define and discuss in detail in Appendix A) are sensi- 
tive to different system scales, and this is especially help- 
ful for spatially-embedded systems like granular packings 
(see Fig.[2|. Our results indicate that the global efficiency 
Ey^ [see Eqs. ([5| and (23)] is a system-level property 
with smallest values along the perimeter of the system 
and largest values in the center. Community structure 
and its associated community label X [see Eqs. ([l5| and 
(27)] and intra-community strength z-score [see Eqs. |3] 
is a meso-scale property and can be used to probe inter- 
mediate structural features. We find that geodesic node 
betweenness B^^ [see Eqs. (IgI) and (25)] can be thought of 



as a one dimensional property in these materials because 



Identifying a Characteristic Size Scale 

An ongoing challenge in the study of granular systems 
is identifying and measuring characteristic size scales 
within granular materials, from the perspective of either 
particles or force chains ^ • Network modularity 
provides a novel means to measure such size scales via the 
identification of community sizes. We find that the op- 
timal value of modularity is a reliable diagnostic for the 
structure of both the binary and weighted networks (see 
Appendix B). To seek characteristic community sizes, we 
also examine community structure as a function of a res- 
olution parameter 7 ^34H36] . The modularity index is 
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FIG. 3. (A,C) Modularity index Q^^ and (B,D) number of 
communities as a function of the resolution parameter (A,B) 
7 and (C,D) the effective fraction of antiferromagnetic edges 
^(7). Error bars indicate the standard deviation over the 17 
experimental runs. 



where A^in is the largest number of negative entries 
Jij (7) for which an A^-node network forms a single com- 
munity and Amax is the smallest number for which the 
network still forms N communities of size 1. 

We examine community structure as a function of the 
transformed resolution parameter ^(7), which we vary 
between and 1. The optimized modularity and the 
number of communities both change gradually for most of 
the ^(7) range (see Fig.[3^,D), although abrupt changes 
are evident for very low and very high values of ^(7). The 
gradual change hints at an interesting size scale, which 
occurs in partitions that contain about 50 — 100 commu- 
nities (with a characteristic size of roughly 5 — 8 parti- 
cles). One possibility is that this size corresponds to the 
width of a shear band, which arise in a variety of mate- 
rials with particulate structure Another possibility 
is that this size corresponds to the 'cutting' length scale 
[43], which is set by a community size at which the 
excess (over-constrained) number of contacts in the bulk 
of a region is equal to the number of contacts around 
the perimeter. If this latter association is correct, then 
the mean number of particles per community would scale 
with the confining pressure. Future experiments can test 
this hypothesis. 



m 

Qw = "^[Wij - -fPij]S{gi,gj) , (1) 

where node i is assigned to community node j is as- 
signed to community gj^ ^{gi^gj) = 1 if gi = gj and it 
equals otherwise, and Pij is the expected weight of the 
edge connecting node i and node j under a specified null 
model. We used the usual Newman-Girvan null model, in 
which the expected strength distribution of the network 
is preserved but ends of edges are rewired uniformly at 
random [ST.'SF. We employed the Louvain locally greedy 
algorithm to optimize modularity [39 , and we varied the 
resolution parameter 7 from 0.001 to 100. Low values of 
7 probe large spatial scales, and high values probe small 
scales. When we increase 7, the number of communities 
increases (as expected), and the modularity decreases. 
See Fig.[3]A,C. 

One can think of the term Jij (7) = Wij — jPij in 
Eq. as a particular choice of interaction strength be- 
tween a pair of spins in a Potts model [34, 37, i40j. We 
exploit this analogy with the Potts model to transform 
the resolution parameter 7 so that it measures the effec- 
tive fraction of antiferromagnetic edges ^(7) in a network 
[41]. We define ^^(7) to be the number of negative ele- 
ments of J(7). The transformed resolution parameter is 
then 



Geography of Community Structures 

Using modularity optimization j36H38l |44] , we find 
that the force-chain network exhibits geographically- 
constrained community structure: groups of particles in 
close spatial proximity are more likely to be a part of the 
same community (i.e., to contact one another with a large 
force) than particles that are farther apart. We examine 
this local neighborhood structure over a variety of size 
scales by varying the resolution parameter 7. We show 
representative results for large spatial scales in Fig. |4]A- 
C. We also note that the communities that we identify in 
granular force networks resemble those in spatial entities 
like states or countries, whose borders are determined 
in part by physical boundaries between neighboring geo- 
graphic domains. 

Importantly, because the optimization of Q is NP-hard 
[45 , one does not expect an optimization algorithm to 
give a global optimum of Q. Instead, we harness numer- 
ous near-degeneracies among good local optima of 
Q by estimating Q 100 times. We find that the these 
100 values of Q for a given run at a given 7 vary by 
approximately 1 x 10~^^, and the similarity in particle 
assignments to communities is approximately 0.98. We 
quantify this using partition similarity [47l |48], which 
ranges from (not similar at all) to 1 (identical). These 
results indicate that the local geographic structures that 
we are identifying in the 2D granular system are robust, 
suggesting the potential for identifying reproducible 2D 
'geographic' regions. 
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To probe the role of each particle in the community 
structure of a force-chain network (see Fig. |4)l)), we use 
the intra- community strength z-score Zi to measure how 
well connected a node is to other particles in its com- 
munity and the participation coefficient Pi to measure 
how the connections emanating from a particle are spread 
among particles in the different communities [49] . 

The intra-community strength z-score is 



Go 

where Sg. denotes the strength (i.e., total edge weight) 
of node z's edges to other nodes in its own community 
the quantity Sg. is the mean of Sg. over all of the nodes 
in and (JSg. is the standard deviation of Sg. in gi. The 
strength of node i is denoted by Si and gives the total 
force of all contacts on particle i. 
The participation coefficient is ^49J 

p. = 1-1: (I) . w 

where Sig is the strength of edges of node i to nodes in 
community g [49 . 

In Fig. [4^-F, we show the intra-community strength 
z-score and participation coefficient for the community 
structure depicted in Fig. |4]A.. Particles that are geo- 
graphically central to a community tend to have higher 
values of Zi and lower values of Pi than particles at the 
geographic periphery of a community. From a physi- 
cal perspective, Zi tends to be highest in particles with 
many force chains passing through them (compare, e.g.. 
Figs. |4)l) and [4]3) and high values of Pi are associated 
with the boundaries between communities (where there 
are few force chains). 

We test whether the observed properties of commu- 
nity structure in our granular systems are related statisti- 
cally to the inter-particle forces that constitute the force- 
chain structure (see Fig. [4)l),G) by examining the re- 
lationship between intra-community strength z-score z^, 
participation coefficient P^, the normalized node strength 
SI = Si/N^ (i.e., the mean force of all edges emanating 
from a node) and the amplitude of the acoustic signal AI 
using the Spearman rank correlation coefficient p, which 
is defined as the Pearson correlation coefficient between 
ranked variables. We use the Spearman coefficient rather 
than the Pearson coefficient due to the non-normal dis- 
tributions of AI values over particles. 

We find that the mean force of all contacts on a particle 
(i.e., S^) is significantly positively correlated with Zi (see 
Fig. [4^) . The mean value of the Spearman rank correla- 
tion coefficient p over experimental runs and resolution- 
parameter values is p ~ 0.71 ±0.02 (where 0.02 gives the 
standard deviation over experimental runs). This strong 
positive correlation indicates that particles at the centers 



of communities are likely to have more or stronger force 
chains running through them. We also find that S^ is 
negatively correlated with P (Fig. |4j). The mean p over 
experimental runs and values of the resolution parameter 
is p ~ —0.05 ± 0.04, where we again take the standard 
deviation only over experimental runs. This negative cor- 
relation indicates that inter-community boundaries occur 
at particles with fewer or weaker force chains. Note addi- 
tionally that a large fraction of the particles have P = 0. 
This is a consequence of the fact that the communities 
are geographically constrained such that the majority of 
particles have contacts only within their own community. 

The relationship between z, P, and S' is expected 
mathematically. For example, if the edges of node i 
all lie within its own community, then S^ and z are re- 
lated linearly according to the following equation: Si = 

X <^Sg. ~^Sg^ , where Sg. is the strength of edges of node i 
to other nodes in its community ^i, and Sg^ is the mean of 
Sg. over all of the nodes in gi. This linear relationship is 
evident for the four communities that we show in Fig.[4]H[. 
Nodes whose connections span more than one community 
(so-called 'boundary nodes', for which the value of P is 
greater than 0) are not so simply related. 

Signal Propagation on Force- Weighted Contact 
Networks 

Previous work in Ref. [22] has shown that the propa- 
gation of acoustic signals is facilitated along strong force 
chains in granular materials, via the increased contact 
area at strong contacts. With this in mind, we test 
whether the geographic community structure of force- 
chain networks is related to signal propagation. As the 
example shown in Fig. [5]A. indicates, we found that 
which we measured over a range of size scales associ- 
ated with resolution-parameter values 7 G [0.001,100], 
is significantly correlated with the signal amplitude AI. 
(For this example run, p ^ 0.57 and the p-value is 
p ^ 2.1 X 10~^^.) The statistical correlation between 
network structure and signal amplitude exists not only in 
the highly heterogeneous signal injection phase, in which 
sound propagates from the driver to nearby particles, but 
also in the more homogeneous scattering phase, in which 
sound reverberates throughout the system. In Fig. |5^, 
we show the results at 7 = 0.1; Figure |5p shows that the 
mean p over runs, 7 values, and time is 0.34 ±0.06. This 
suggests that similar dynamic principles underlie sound 
propagation in both injection and scattering phases. We 
discuss the dynamics within the two phases in more detail 
in the next section. 

In quantifying the relationship between the intra- 
community strength z-score z (a property of the algo- 
rithmically computed community structure) and the sig- 
nal propagation amplitude A/, we note that the corre- 
lation between these two variables is decreases as 7 in- 
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FIG. 4. [Color online] The geographic sizes of communities tend to decrease as the resolution parameter 7 is increased from 
(A) low (7 = 0.1 to (C) high 7 = 1 values. We color particles according to their community label. One can examine community 
structure of the (D) force-chain network using geographic location, (E) intra-community strength z-score Zi, (F) participation 
coefficient Pi, and (G) mean force per particle Si. Example scatterplots for a single experimental run showing that mean force 
S'i on particle i is (H) positively correlated with the intra-community strength z-score (the Spearman correlation coefficient 
for this run is p ^ 0.96 and the p-value is p ^ 1.9 x 10~^^^) and (I) negatively correlated with the participation coefficient 
(p ~ —0.18 and p ~ 1.8 x 10~^). In panel (H), nodes assigned to communities 1 through 4 and whose participation coefficients 
are equal to are displayed using different colored markers. Nodes in any of the 4 communities whose participation coefficients 
are greater than (so-called 'boundary nodes') are displayed using purple markers. The mean correlations for z and P over 
experimental runs and values of the resolution parameter are p ~ 0.87 di 0.08 and p ~ —0.16 di 0.05, respectively. We show the 
results for one experimental run (#2) in this figure, and the results for the other runs are similar. 



creases (see Fig. ^p). The strength of the relationship 
between network structure and signal propagation for 
small 7 suggests that partitions with a few large com- 
munities achieve better estimates of the propagation be- 
havior. Indeed, as demonstrated in Fig. |5jl), when the 
network forms a single community (for very low 7 val- 
ues), the correlation between z and AI is similar to that 
between the mean force per particle {S^) and AI. 

The retention of a correlation between z and AI for 
larger values of 7, for which the network is partitioned 



into more (and smaller) communities, stems from the 
strong correlation between intra-community strength z- 
score and the mean force (normalized strength) of a 
particle (see Fig|4]H[), the latter of which is a particle- 
scale measurement and is independent of spatial resolu- 
tion. The relationship between the meso-scale (commu- 
nity structure) and particle-scale (mean force on a parti- 
cle, which is equal to a node's normalized strength) net- 
work properties stems from the physical embedding of 
the granular system in R^. A particle that is located 
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FIG. 5. [Color online] (A) An example scatterplot between 
logarithms of the intra-community strength z-score [log]^Q(^ + 
5)] and the amplitude of the acoustic signal [log^Q(A/] at 
t — 1 for a single experimental run [j — 2) and resolution- 
parameter value (7 = 0.1). The constant 5 was added to z to 
ensure that all values were positive prior to taking the loga- 
rithm. The Spearman rank correlation coefficient is p ~ 0.57 
and the p- value is p ^ 2.1 x 10~^^. (B) Correlation between 
logio(zH-5) and logio(AI) for all 17 runs as a function of 
time: t < 40 is the acoustic signal injection phase, and t > 40 
is the acoustic signal scattering phase. In the bottom part 
of panel (B), we show a trace of the voltage 1/ of a piezo 
particle; the injected signal has an amplitude of 0.5V. (C) 
Correlation between logio(z + 5) and logio(AI) as a function 
of 7, where 7 G [0,2] is increased in increments of 0.1. We 
averaged the correlation over all 80 time points at which A/ 
was measured, and the mean p over runs 7 values, and time 
was 0.34 ± 0.06. Box plots show the variability over experi- 
mental runs. (D) Correlation between logio(AI) and a vari- 
ety of network diagnostics: intra-community strength z-score 
z — ^(7) for 7 = 0.001, 7 = 10, and 7 = 100; and weighted 
(white background, left of the figure) and unweighted (light 
gray background, middle of the figure) versions of global ef- 
ficiency [Eyj{i) and E{i)\^ geodesic node betweenness [Byj(i) 
and 5(z)], and clustering coefficient [Cw{i) and C{i)\. For 
completeness, we also show results for the mean force S' (left 
of the figure), which was the variable previously reported to 
be correlated with A/ l22l. 



geographically inside of a community has all of its con- 
nections to other particles in its community because it is 
constrained to connect only to its geographic neighbors 
(i.e., there can be no long-range contacts). This is unlike 
most investigated real- world networks [36l [37] , in which 
communities tend to be highly interconnected and most 
nodes have at least some connections to nodes in other 
communities. 

To assess whether community structure is unique in its 
ability to predict signal amplitude, we also examine other 



weighted network diagnostics that are sensitive to differ- 
ent system dimensionalities (see Fig. [2| . Our results sug- 
gest that community structure (see Fig. |2^) is a better 
predictor of signal propagation than system-scale (e.g., 
global efficiency; see Fig. [2j^), curve-scale (e.g., geodesic 
node betweenness; see Fig. [2]C), and particle-scale (e.g., 
clustering coefficient; see Fig. |2|l)) network diagnostics. 
See Appendix A for mathematical definitions and intu- 
itive descriptions. In particular, clustering coefficient and 
A/ are not strongly correlated, so triangles of contacts 
do not appear to be important for signal propagation. 
See the comparison in Fig. [5p). 



Phase Sensitivity of Network Diagnostics 

Although the correlation between z and signal ampli- 
tude is strong in both injection and scattering phases for 
small 7 (i.e., large community size), it is higher in the 
scattering phase than in the injection phase when aver- 
aged over all resolutions (7 G [0.001, 100]; see Fig.|6^,E) 
We do not observe such sensitivity to phase for the mean 
force per particle (see Fig.|6|^,D). Interestingly, the signal 
propagation during the injection phase is more strongly 
correlated with the global efficiency than it is during 
the scattering phase (see Fig.|6]C,F), suggesting that the 
acoustic signal propagates over the shortest weighted 
paths during the injection phase. These results illus- 
trate insights from network analysis that one cannot ob- 
tain from particle-scale measurements: signal propaga- 
tion during injection is well characterized by shortest 
paths that span the system, whereas it is characterized 
by local neighborhood structure during scattering. An 
interesting question is whether the amplitude of the in- 
jected signal affects the size of the geographic neighbor- 
hood through it propagates. 

We also examine the sensitivity of the relationship be- 
tween z and A/ to the injection and scattering phases 
as a function of the resolution parameter (see Fig. [7]A). 
The correlation between z and signal amplitude is con- 
sistently higher in the scattering phase than in the injec- 
tion phase throughout 7 G [0.001, 100]. Furthermore, the 
largest difference in the Spearman correlation between z 
and A/ for the scattering versus injection phases occurs 
for partitions with approximately 50 — 100 communities, 
corresponding to community sizes of roughly 5 — 8 par- 
ticles (see Fig. [7^). This is similar to the size scale that 
we identified previously when using the transformed res- 
olution parameter. 



DISCUSSION 

A networks perspective provides a useful framework 
in which to study the material and dynamic properties 
of granular materials. Network diagnostics vary in their 
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FIG. 6. [Color online] Spearman correlations between signal 
amplitude AI and (A) the mean force per particle S\ (B) 
the intra- community strength z-score^; for 7 = 2, and (C) the 
global efficiency . Data is for all experimental runs and for 
all times, including both injection (left) and scattering (right) 
phases. Box plots of the Spearman correlations between signal 
amplitude AI and the three diagnostics shown in (A-C): (D) 
mean force per particle (ps')^ (E) intra-community strength 
z-score z (pz), and (F) global efficiency (pE^)- We have aver- 
aged the correlations that we show in the box plots over the 
injection (left) and scattering (right) phases. For (E), note 
that we also average the correlations over the resolution pa- 
rameter 7. Using MATLAB notation, the precise values of 7 
that we considered are [0.001 : 0.001 : 0.009, 0.01 : 0.1 : 1, 2 : 
3, 4 : 0.1 : 20, 30 : 10 : 100, 200 : 100 : 1000]. The reported 
p- values indicate the results of 2-sample t-tests. 



sensitivity to scales of the granular system: the particle- 
scale can be probed with a clustering coefficient, the 
curve-scale can be probed with geodesic node between- 
ness, the domain-scale can be probed with community 
structure, and the system-scale can be probed with global 
efficiency. Moreover, one can identify potentially inter- 
esting length/size scales in the system using meso-scale 
network features such as community structure. As we 
show in Appendix C, one can also obtain physical in- 
sights into the geographic organization of the material 
by comparing the features of the actual networks to a 



0.5 
0.45 

0.4 
I 0.35 
I 0.3 

o 

^ 0.25 
0.2 




50 

resolution parameter, y 




100 200 300 

number of communities 



FIG. 7. [Color online] (A) Spearman correlations between 
signal amplitude AI and the intra-community strength z- 
score z averaged over the injection (black) and scattering 
(red) phases as a function of the resolution parameter 7. (B) 
Difference between the correlation in the scattering (pi) and 
injection (pi) phases as a function of the number of communi- 
ties in the partition. We find the greatest sensitivity to phase 
for partitions with approximately 50 — 100 communities (i.e., 
for communities containing roughly 5 — 8 particles), which 
corresponds to the size scale that identified as potentially in- 
teresting using the transformed resolution parameter. 



null model consisting of an ensemble of random geomet- 
ric graphs. 

The dynamics of signal propagation on a network are 
best characterized by weighted diagnostics derived from 
the granular force-chain network, suggesting that the 
topology of the underlying (unweighted) contact network 
alone is not sufficient to explain signal propagation. In 
other words, one must also consider network geometry. 
This result underscores the important relationship be- 
tween signal propagation and force-chain organization 
(see Fig.jSp). Similar phenomena are likely relevant for a 
variety of energy transport problems (e.g., in sound, heat, 
and electricity) in a broad class of amorphous materials. 

The algorithmic detection of communities is particu- 
larly useful in quantifying the effect of meso-scale net- 
work structure on signal propagation. We find that com- 
munity structure is a better predictor of signal ampli- 
tude over the range of propagation phases than system- 
scale, curve-scale, or particle-scale measurements. Fur- 
thermore, by contrasting signal behavior during injection 
and scattering phases, we are able differentiate the sensi- 
tivities of system-wide and meso-scale network structure 
to sound propagation. Community structure seems to be 
a better predictor of signal propagation in the scattering 
phase than in the injection phase, suggesting that the 
sound scatters in local geographic neighborhoods. How- 
ever, global efficiency predicts signal propagation bet- 
ter in the injection phase than in the scattering phase, 
suggesting a more system-wide dynamic distribution of 
sound. 

Studying community structure allows one to investi- 
gate the meso-scale architecture of the granular pack- 
ing. In addition to the intra-community strength z- 
score predicting particle sound amplitude, the identifi- 
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cation of geographic communities provides a quantifiable 
size scale, which may be useful for seeking the diverg- 
ing length scale that is expected at the jamming transi- 
tion . The meso-scale nature of the sound propaga- 
tion might be related to other meso-scale phenomena in 
granular physics, such as the spatial eigenmodes for soft 
(low-energy) modes, which are observed in simulations to 
take the form of localized swirls [50l [51] . Our approach 
also provides a framework to relate ID structures (force 
chains) to 2D structures (geographic domains). It there- 
fore might also prove useful in other settings — such as 
in the study of crystalline solids, where domain structure 
is critical to system function. 

The presence of correlated regions such as geographic 
communities in a granular material is reminiscent of 
shear transformation zones (STZs [52 ), in which local- 
ized regions throughout a sheared material have a higher 
propensity to deform under shear. Importantly, how- 
ever, the community structure that we compute spans 
the system, whereas STZs are relatively small structures 
dispersed throughout system. Also of interest is a com- 
parison with the results of Ref. |53], which illustrated 
that vibrational modes can identify soft spots in sheared 
systems. 

In conclusion, using network analysis to study granu- 
lar materials can be extremely useful, as it can help char- 
acterize particle, curve, domain, and system-scale prop- 
erties of such materials. In particular, the algorithmic 
detection of communities provides a means to identify 
potentially interesting characteristic size scales in such 
systems. When combined with time-resolved acoustic 
measurements [22,, such a networks perspective can il- 
luminate the meso-scale structures within which sound 
travels preferentially. We found that particles that are 
well connected to their community have larger-amplitude 
signals passing through them. Our results also suggest 
that signals scatter in local geographic neighborhoods 
but propagate more systemically during signal injection. 
Investigation of both weighted and unweighted networks 
demonstrates that a weighted network is a better predic- 
tor of sound propagation, suggesting that the force-chain 
structure of the granular material is an important compo- 
nent in sound propagation. Our results demonstrate that 
one cannot examine only system-scale or local-scale net- 
work features to understand how sound travels through 
a granular material. Importantly, one achieves a better 
description of sound propagation when one includes how 
the particles relate to their neighbors in a network. 
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APPENDIX A: DEFINITIONS OF NETWORK 
DIAGNOSTICS 

Diagnostics Applied to Unweighted Contact 
Networks 

To characterize structure of the binary (contact) net- 
works, we examined 21 diagnostics: number of nodes, 
number of edges, global efficiency [54 , geodesic node be- 
tweenness centrality ^1 , random- walk node between- 
ness [57], geodesic edge betweenness [58 , eigenvector cen- 
trality [59^, closeness centrality [60 , subgraph centrality 
[61 , communicability [62], clustering coefficient [63 , lo- 
cal efficiency [54] , modularity optimized using two differ- 
ent algorithms [37^ '39', '4T| , hierarchy [64^, synchronizabil- 
ity [65j, degree assortativity [66 , robustness to targeted 
and random attacks [67], the Rent exponent [68], and 
mean connection distance [69^. 

In our descriptions below, we give for each diagnostic 
(i) a mathematical definition, (ii) an intuitive description 
of the term, and (iii) a comment on its possible physical 
significance for the granular system that we study. We 
also computed node-specific values for the following di- 
agnostics: geodesic betweenness, global efficiency, and 
clustering coefficient. (See the discussions below.) 

1. Number of nodes N: (i) The diagnostic TV is defined 
as the number of nodes in a network, (ii) It is used 
as a measure of the size of a system, (iii) In this 
study, N is the number of particles in the system. 
It provides a consistent but dynamically uninter- 
esting characterization of the network because it is 
identical at all points in time. 

2. Number of edges D: (i) The diagnostic D is defined 
diS D = "^.j Aij , where A is an unweighted (binary) 
network with components A^j. Nodes are particles, 
and an edge exists between particles i and j (i.e, 
Aij = 1) if and only if particles i and j are in con- 
tact with each other (otherwise, Aij = 0). (ii) The 
quantity D is simply the total number of edges in 
the system, (iii) The number of edges D is related 
to the mean contact number, which is denoted by 
z in the granular-materials community. The mean 
contact number of the system is equal to z = 
The diagnostic D provides a consistent but uninter- 
esting characterization of the network because the 
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number of contacts scales with pressure [2] (which 
is the same for ah experimental runs). 

3. Global efficiency E [54 : (i) Let dij be the shortest 
(geodesic) number of steps necessary to get from 
node i to node j. The global efficiency is then de- 
fined as 

(ii) Global efficiency can be interpreted as a mea- 
sure of how well a signal is transmitted through 
a network, (iii) One can expect the global effi- 
ciency to be small in 2D granular packings because 
particles that are not geographically close to one 
another are separated by multiple contacts (edges) 
and therefore by a long path length (low efficiency). 
As one can see in Table \\T\\ this is indeed the case. 

4. Geodesic node betweenness B [55 : (i) Geodesic 
node betweenness is defined for the i^^ node in a 
network Q as 

B,= Z *f^, (6) 

where all three nodes (j, m, and i) must be differ- 
ent from each other, V^j,m is the number of geodesic 
paths between nodes j and m, and il^j^mii) is the 
number of geodesic paths between j and m that 
pass through node i. The geodesic betweenness of 
an entire network B is defined as the mean of Bi 
over all nodes i in the network, (ii) Geodesic be- 
tweenness can be interpreted as a measure of traffic 
flow on a network, (iii) One might expect the ma- 
jority of geodesic paths that link any node of the 
packing to any other node to pass through the mid- 
dle of the system. Indeed, we find that the largest 
values of betweenness occur in the center of the sys- 
tem and the smallest values along the edges of the 
packing. 

5. Random- walk node betweenness B^w [Si]- (i) For 
an adjacency matrix A and diagonal matrix D, let 

= At • D^^ be the matrix M with the row 
and column t removed (and A^ and are defined 
analogously). The probability that a walk starts 
at s, takes n steps, and ends up at some node i 
(which cannot be t because t has been removed) 
is given by element is of MJ^; denote this element 
by [MJ^]i5. Walks end up at v and w with prob- 
abilities [MJ^]^5 and [MJ^J^^. Fractions l/k^ and 
1/kw of these walks subsequently pass along the 
edge {v^w) in one direction or the other, assuming 
that such an edge exists. (Note that ky is the de- 
gree of V and k^ is the degree oiw.) Summing over 
all n shows that the mean number of times that a 



walk of any length traverses the edge from v to w 
is k~^[{l — yit)~^]vs' The random- walk between- 
ness of a node is the mean of this quantity over all 
edges emanating from that node, and the random- 
walk betweenness of the entire network is the mean 
of the random- walk betweenness of all nodes in the 
network, (ii) Random-walk betweenness can be in- 
terpreted as a measure of information flow or signal 
flow in a network, (iii) Similar to geodesic node be- 
tweenness, one might expect the random- walk node 
betweenness to be highest in the center of the sys- 
tem and lowest on the edges of the system. This is 
indeed the case. 

6. Geodesic edge betweenness B^ [W: (i) Inspired by 
Freeman's geodesic node betweenness, the geodesic 
edge betweenness of an edge is defined as the num- 
ber of shortest paths between pairs of nodes that 
run along it. For the edge connecting nodes j and 
m, the geodesic edge betweenness is given by 

Be ( j, m) = ^ ( j, m) , (7) 

where i/ji^kU^'^^) is the number of shortest paths 
between i and k that pass through the edge con- 
necting nodes j and m. (ii) One can interpret edge 
betweenness as a measure of the influence of an 
edge on traffic flow through a network, (iii) In a 2D 
granular packing, edge betweenness might indicate 
the influence of a contact on a hypothetical flow 
through the network. In our system, we find that 
edge betweenness is largest in the center of the sys- 
tem and smallest on the edges of the system. This 
is consistent with the results for the geodesic node 
betweenness. 

7. Eigenvector centrality Ce [70]: (i) The eigenvector 
centrality Ce{i) of node i is proportional to the sum 
of the centralities of the nodes connected to it: 

Ce{^) = \ E CeU) = \^AMj). (8) 

jeMii) j 

where M{i) is the set of nodes that are neighbors 
of i (i.e., which are connected to i directly via an 
edge) and A is the largest eigenvalue of A. From 
equation ([8|, one can deduce that Ce{i) is the i*^ 
component of the leading eigenvector (each entry of 
which is positive by the Perron-Frobenius theorem 
[1 ) of the adjacency matrix, (ii) Eigenvector cen- 
trality can be used to measure the importance of 
a node in a network based on its direct connection 
to important nodes, (iii) In a 2D granular pack- 
ing, one would expect eigenvector centrality to be 
large for a particle that has many contacts or for 
a particle whose immediate neighbors have many 
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contacts. Indeed, we find that eigenvector central- 
ity is iiigiiest in a local region of the system in which 
high-degree nodes are most concentrated. 

8. Closeness centrality Cc [60 : (i) We use a version 
of closeness centrality that is appropriate for both 
connected and disconnected graphs [71 . It is de- 
fined as 



Cc{i) 



jev/i 



(9) 



where il^g{i^j) is the geodesic distance between 
nodes i and j (i.e., the length of the shortest path 
connecting i and j) and the notation V/i indicates 
that V is the connected network component reach- 
able from i and does not include i. (ii) Closeness 
centrality can be used as a measure of the impor- 
tance of a node in a network, (iii) For 2D granular 
systems, one might expect closeness centrality to 
be small given the lattice-like topology of a contact 
network. However, as shown in Table |TITj closeness 
values are somewhat larger than those for random 
geometric graphs (RGGs; see the discussion in Ap- 
pendix [] 

9. Subgraph centrality Cs [HP: (i) We first note that 
the number of closed walks of length k starting and 
ending at node i is given by the k^^ local spectral 
moment which is defined as the i^^ diagonal 

entry of the k^^ power of the adjacency matrix A: 



(10) 



The subgraph centrality of node i is then defined 
as 

^^.w = E^- (11) 



k=0 



k\ 



(ii) Subgraph centrality characterizes the participa- 
tion of each node in all subgraphs in a network, (iii) 
For 2D granular systems, one might expect sub- 
graph centrality to be small because nodes partici- 
pate in few subgraphs other than their own, as their 
connectivity is strongly constrained to their local 
geographic neighborhood. Indeed, as indicated in 
Table |III[ the subgraph centrality has a value that 
is less than 1/3 that of the value that we computed 
for a corresponding ensemble of RGGs (see our later 
discussion). 

10. Communicability Co [62]: Because of the direct 
relationship between the powers of the adjacency 
matrix A and the number of walks in a network, 
one can define the communicability between nodes 
i and j as 



Coij 



E 

/c=0 



k\ 



(12) 



The communicability Co of a network is then the 
mean of the communicabilities of each pair of (non- 
identical) nodes, (ii) Communicability was devel- 
oped to measure the ease of communication or 
transmission in terms of passage between different 
nodes in a network, and it is specifically based on 
walks rather than paths [62 . (iii) For 2D granular 
systems, one might expect the mean communica- 
bility to be small because the geographic nature of 
the contacts creates a lattice-like topology. Indeed, 



as indicated in Table III, its value is less than 1/4 



than the value that we computed for a correspond- 
ing ensemble of RGGs. 

11. Clustering coefficient C [63]: (i) The diagnostic C 
is defined by supposing that a node i has ki neigh- 
bors, so a maximum oi ki{ki — 1) /2 edges can exist 
between these neighbors. The local clustering co- 
efficient Ci is the fraction of these possible edges 
that actually exist: 



Ci 



ki{ki 1) 



(13) 



The clustering coefficient C of an entire network is 
then defined as the mean of Ci over all nodes i. (ii) 
The clustering coefficient C can be interpreted as a 
measure of local clustering properties in a network, 
(iii) One can expect C to be large in 2D granular 
packings because particles that are geographically 
close to one another are also near each other in a 
network. This ought to yield a large number of 
connected triples and hence a high value of C As 
shown in Table [III] we do indeed observe reasonably 
large values [72 for clustering coefficients in the 17 
experimental runs (the mean value over all runs is 
C ~ 0.26), but interestingly the mean value of C in 
the corresponding RGG ensemble is twice as high. 

12. Local efficiency Ei [54]: (i) The the local efficiency 
of node i is defined as 



Eiii) = 



(14) 



where Gi is the subgraph consisting of all nodes 
connected to node i along with all of their edges 
between each other, and dj^^ is the minimum path 
length between nodes j and k in this subgraph. The 
local efficiency Ei is the mean value of Ei{i) over 
all nodes i. (ii) Local efficiency Ei can be inter- 
preted as a measure how well a signal is transmitted 
through a subgraph, (iii) One might expect local ef- 
ficiency to be large in 2D granular packings because 
particles that are very close to each other geograph- 
ically lie in one another's subgraphs. However, as 
we show in Table [lll| we obtain values that are only 
about half of those for corresponding RGGs. The 
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mean granular- network value of 0.33 is comparable 
in value to some communication networks ^54^ . 

13. Modularity index Q [36^-^ : (i) Networks can be 
partitioned into communities (or modules) in which 
nodes inside the same community are more densely 
connected to each other than they are to nodes in 
other communities. The modularity of a network 
partition is defined as 



Q=—y 



A. 



'2D 



(15) 



where ki is the degree of node D is the total num- 
ber of edges in the network, 5ij is the Kronecker 
delta, and gi is the community to which node i 
has been assigned. With the standard null model 
Pij = kikj /{2D)^ equation (15) is sometimes called 
' Newman- Girvan modularity.' One uses one of nu- 
merous possible computational heuristics to maxi- 
mize Q in the space of all network partitions, and 
one can then report the maximum value obtained 
for Q. However, it is important to note that the 
optimization of Q is NP-hard [45], so one cannot 
expect the output of an optimization algorithm to 
be a globally optimal partition of a network. In this 
light, we use two different computational heuristics 
to optimize Q: Newman's spectral algorithm [4? 
(which yields a modularity value that we denote 
Qs) and the Louvain locally greedy method [39] 
(yielding a modularity value that we denote Ql)- 
(ii) The optimal value of Q is a measure of how well 
a network can be partitioned into cohesive commu- 
nities, (iii) In a 2D granular system, one might ex- 
pect communities to be localized geographically be- 
cause connectivity between nodes in potential com- 
munities is constrained geographically. Indeed, as 
shown in Table |III[ the values of Qs and Q l are 
both extremely high [36l |37] . 



14. Hierarchy h [64]: (i) A sense of hierarchical struc- 
ture in a network can be characterized by the coeffi- 
cient /i, which is used to quantify a putative power- 
law relationship between clustering coefficient Ci 
and the degree ki of all nodes in the network [64] : 



Ci 



(16) 



(ii) Networks in which clustering coefficient has a 
power-law scaling with degree possess a hierarchy 
in which hubs (i.e., high-degree nodes) tend to have 
low local clustering and low-degree nodes tend to 
have high local clustering. The parameter h gives a 



precise scaling of such effects when the relation ( 16 ) 



holds, and it can perhaps indicate a looser sense of 
hierarchy in more general situations, (iii) It is not 
clear a priori whether 2D granular packings have 



some hierarchical characteristics, though the au- 
thors of Ref. [64 have suggested that geographic 
networks are not very hierarchical. Our calcula- 



tions (see Table III) provide some support for this 
claim, as we observe significant scaling between Ci 
and ki and obtain h ^ 0.76. It is noteworthy, how- 
ever, that this value is roughly three times what 
one obtains in a corresponding RGG ensemble (see 
Table [Hi). 



15. Synchronizability s 
is defined as 



(i) The synchronizability 



(17) 



where A2 is the second smallest eigenvalue of the 
Laplacian C of the adjacency matrix and A at is the 
largest eigenvalue of C [65 . (ii) The synchronizabil- 
ity of a network characterizes structural properties 
of a network that hypothetically enable it to syn- 
chronize rapidly, (iii) One might expect that the 
synchronizability of the contact network in a 2D 
granular packing is small due to the lattice-like na- 
ture of the network topology. Indeed, as shown in 
Table [Hlj the value for s for our system is tiny. 



16. Degree assortativity a [66 : (i) The degree assor- 
tativity of a network (which is often called simply 
'assortativity') is defined as 



(18) 



where ji and ki are the degrees of the nodes at 

the two ends of the i*^ edge (i e {1 E} [44]. 

(ii) Degree assortativity measures the preference of 
a node to connect to other nodes of similar degree 
(leading to an assortative network, for which a > 0) 
or to nodes of very different degree (leading to a dis- 
assortative network, for which a < 0). (iii) It is not 
clear a priori whether 2D granular packings should 
display degree assortativity. Our calculations indi- 
cate that the degrees exhibit some mild positive as- 
sortativity (a ~ 0.14), but the corresponding RGG 
ensembles have a significantly higher positive as- 
sortativity ofa^ 0.56. (See Table [Hi]) 



17. Robustness R [67|: (i) One can define robustness 
in terms of different types of attacks on a net- 
work. In the most commonly studied type of tar- 
geted attack, nodes are removed (one by one) in 
descending order of their degree; in a random at- 
tack, nodes are removed in random order. Each 
time a node is removed from a network, we re- 
calculate the size S (i.e., number of nodes) of the 
largest connected component. One can examine ro- 
bustness by plotting S versus the number of nodes 
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removed n [73H75] . One can then define a robust- 
ness parameter R as the area under the curve in 
the plot of S = S{n). More robust networks retain 
a larger connected component even when several 
nodes have been removed; this is represented by 
a larger area under the curve and hence by larger 
values of R. (ii) Robustness indicates network re- 
silience to either targeted {Rt) or random (Rr) at- 
tacks, (iii) Robustness tends to be most interesting 
for networks with highly heterogeneous degree dis- 
tributions, so it might not be very insightful for 
2D granular packings. We note, however, that we 
find values of Rt and Rr for our granular networks 
that are more than 3 times as large as those for the 
corresponding RGG ensemble. (See Table [lll[) 



18. Rent exponent p [76 : (i) Rent's rule, which was 
first discovered in relation to computer chip design, 
defines a scaling relationship between the number 
of external signal connections (edges) e to a block 
of logic and the number of connected nodes n in 
the block [76]: 



e ^ w 



(19) 



where p G [0, 1] is the Rent exponent, (ii) The 
Rent exponent measures the efficiency of the phys- 
ical embedding of a topological structure, (iii) Due 
to the physical constraints of a 2D granular pack- 
ing, we expect to observe Rentian scaling with a 
relatively low Rent exponent (similar to that of a 
lattice). The theoretically expected minimum of 
the Rent's exponent for a 2-dimensional physical 
lattice is pt = 1-1/De [77H79] where De'is the Eu- 
clidean dimension of the space (e.g., 2), and there- 
fore Pt = 0.5, which is consistent with empirical 
results on memory circuits [80 . However, the ex- 
pected value of the Rent exponent might further 
depend on the type of physical lattice under study 
(e.g., a rectangular or hexagonal lattice). 

19. Mean connection distance mcd: (i) An edge's es- 
timated connection distance Lij is defined as the 
Euclidean distance between the centroids of par- 
ticles i and j. (ii) The mean connection distance 
mcd is defined as the mean of all Lij values in the 
network, (iii) The mean connection distance in a 
2D granular packing is related to the number of 
particles, the area of the system, and particle size. 



Diagnostics Applied to Force- Weighted Contact 
Networks 



To characterize the structure of the force-weighted 
networks, we used 8 diagnostics: normalized strength 



[81] [82], diversity [83], path length [51], geodesic node be- 
tweenness [55[ [82] , geodesic edge betweenness [58^ , clus- 
tering coefficient [82], transitivity [85], and optimized 
modularity [38 . 

1. Normalized Strength 5" [8T1[82: (i) The strength of 
node i is given by the column sum of the weighted 
adjacency matrix: 



(20) 



and the strength of an entire weighted network 
W is the mean of Si over all i. The normalized 
strength S' is 



q/ _ 



(21) 



where N is the total number of nodes, (ii) Strength 
is a measure of how strong the connections are in 
a network, (iii) In the present context, normalized 
strength provides a measure of the mean contact 
forces between particles, and we therefore expect 
this diagnostic to be correlated with sound propa- 
gation. Indeed, we observe this in our calculations. 

2. Diversity V [83]: (i) The diversity of node i is de- 
fined as the variance of the edge weights for the set 
of all edges connected that are connected to it. It 
is given by 



1/2 



(22) 



The diversity of an entire weighted network is the 
mean of Vi over all i. (ii) Diversity is a measure of 
the variance of connectivity strengths in a network, 
(iii) In the present context, diversity is a measure 
of the variance of contact forces between particles, 
and it has a high positive correlation with normal- 
ized strength. 

3. Global efficiency [54]: (i) Let dfj = m8ix{Wij) - 
Wij be the weighted shortest path between nodes 
i and j. The global efficiency of node i is then 
defined as 



N 



—y-- 



(23) 



The global efficiency E^, is the mean value of Euj{i) 
over all nodes i. (ii) One can interpret global ef- 
ficiency as a measure of how efficiently a signal is 
transmitted through a network, (iii) We expect the 
global efficiency of the force-weighted contact net- 
work to be large in the center of the packing and 
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small on the edges of the packing because parti- 
cles that are not geographically close to each other 
do not exert forces on one another. Indeed, this is 
what we observe. 

4. Clustering coefficient C^j [22]: (i) One can define a 
weighted clustering coefficient Cw{i) of node i using 
the formula 

1) ^ 



Si {hi 



AijAikAjk , (24) 



where Si is node i's strength, ki is its degree, W is 
the weighted adjacency matrix, and A is the under- 
lying binary adjacency matrix, (ii) The weighted 
clustering coefficient Cw{i) measures the strength 
of local connectivity, (iii) It is constrained by the 
underlying contact network structure, so we expect 
it to have a high positive correlation with the bi- 
nary clustering coefficient C{i). Indeed, the Pear- 
son correlation coefficient between the binary and 
weighted clustering coefficients over the experimen- 
tal runs is r ~ 0.94 (with a p- value of p ~ 2 x 10~^) . 
Both diagnostics tend to attain their highest values 
on the edges of the packing, where nodes' immedi- 
ate neighbors are most likely to also be connected 
to one another. 

5. Geodesic node betweenness By^ [86 : (i) Geodesic 
betweenness is defined for the i^^ node in a network 
Q as 



(25) 



where all three nodes (j, m, and i) must be dif- 
ferent from each other, tl^j^rn denotes the number 
of geodesic weighted paths between nodes j and 
m, and V^j,m(0 denotes the number of geodesic 
weighted paths between j and m that pass through 
node i. (As with weighted global efficiency, the 
weighted shortest path between nodes i and j is 
dfj.) The weighted geodesic betweenness of an en- 
tire network is defined as the mean of B^^i) 
over all nodes i. (ii) One can interpret weighted 
geodesic betweenness as a measure of traffic flow on 
a network, (iii) We expect betweenness in weighted 
networks to correlate positively with strength, just 
as betweenness in binary networks correlates posi- 
tively with degree [52- In a 2D granular packing, 
we expect particles in the center of the system to 
have high values of weighted betweenness because 
more paths must pass through them to connect the 
edges of the system. We indeed find this to be the 
case. 

6. Geodesic edge betweenness B^^ [66l |86]: (i) We 
define geodesic edge betweenness on weighted net- 
works using the number of shortest weighted paths 



between pairs of nodes that run along it. (We again 
determine the path distance between nodes i and j 
using dfj.) For the edge connecting nodes j and m, 
the weighted geodesic edge betweenness is therefore 



Bew ( j, m) = ^ ^i,/c ( j, rn) , 



(26) 



where 2pi^f^{j,m) is the number of shortest paths 
between nodes i and k that pass through the edge 
connecting nodes j and m. (ii) Weighted edge be- 
tweenness indicates the influence of an edge on traf- 
fic flow through a network, (iii) In a 2D granular 
packing, the edge betweenness should give an indi- 
cation of the influence of a contact on a hypothet- 
ical flow through the network. We find that edge 
betweenness is largest in the center of the system 
because more paths must pass through these edges 
to connect all pairs of particles. 

7. Modularity index [SF-'SS^: The weighted mod- 
ularity of a network partition is 



1 

2W 



E 



2W 



-'gi,gj ' 



(27) 



where Si is node z's strength, W is the total 
strength of the edges in a network, Wij is an el- 
ement of the weighted adjacency matrix, Sij is the 
Kronecker delta, and gi is the label of the com- 
munity to which node i has been assigned. As with 
unweighted networks, one uses some computational 
heuristic to find a partition that maximizes Q. As 
with the binary networks, we have used the Louvain 
locally greedy optimization method [39]. (ii) The 
maximum value of Q is a measure of how well a net- 
work can be partitioned into cohesive communities, 
(iii) In a 2D granular packing, in which forces be- 
tween particles are represented as edge weights, we 
expect communities to be localized in space because 
the forces between nodes in potential communities 
are constrained geographically. Indeed, as shown 
in Fig. |2j this is indeed the case. 



Transitivity T [8 
T{i) of node i is 

m 



(i) The weighted transitivity 



kiiki + 1) 



1/3 



(28) 



where we have normahzed weights by the maximum 
edge weight in the matrix: 



max ( 



(29) 



The transitivity T of the entire network is the mean 
of T{i) over all nodes i. (ii) Weighted transitivity 
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is a generalization of the local clustering coefficient 
in unweighted networks in which one computes the 
sum of the weights of edges in a network's triangles 
instead of computing simply the number of trian- 
gles, (iii) We expect weighted transitivity to be 
similar to the weighted clustering coefficient Cw 
Indeed, we find that the two variables are highly 
correlated over experimental runs (with a Pearson 
correlation coefficient of r ~ 0.99 and a p- value of 
p^lx 10-^^). 



We implemented all computational and simple statis- 
tical operations (such as t-tests and correlations) using 
MATLAB® (2009a, The MathWorks Inc., Natick, MA). 
We estimated network diagnostics using a combination of 
in-house software, the Brain Connectivity Toolbox [87] , 
the MATLAB Boost Graph Library, and the fast unfold- 
ing community detection code for the Louvain optimiza- 
tion of modularity [39] from Peter Mucha [88 . 
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APPENDIX B: RELIABILITY OF NETWORK 
STRUCTURE 

We quantify the reliability of each diagnostic by cal- 
culating the coefficient of variation (a normalized mea- 
sure of dispersion) over the 17 experimental runs: CV = 
cr/|/i|, where a is the standard deviation and /i is the 
mean. Values of CV < 0.2 are commonly considered to 
be acceptable, as they indicate that a variable is reliable 
[iSllHal [^. See Table [l| for CV values for all binary 
network diagnostics and Fig. [9^ for a corresponding bar 
graph. Interestingly, reliable diagnostics are dispersed 
among the quantities we considered rather than focused 
on sets of related diagnostics. 

For the force- weighted granular networks, we find lower 
reliability (i.e., higher values of CV) for the diagnostics 
than for the (binary) contact networks. Compare Fig. [8^ 
to Fig. |9^ and Table [ll] to Table |l|. It is possible that 
the reliability is lower in the weighted networks because 
a large ensemble of force-chain networks are consistent 
with a given packing [9T1. Therefore, for each binary 
network, we are sampling one weighted network out of 
the many possible force-chain networks that could arise 
from the underlying contacts. Based on this degeneracy, 
we might expect that network diagnostics that depend 
on the force topology might be less consistent across ex- 
periments than those based on contacts alone. As with 
the unweighted networks, the strongly reliable weighted- 
network diagnostics are dispersed among the 8 diagnos- 
tics rather than focused on sets of related quantities. 



FIG. 8. [Color online] (A) Relationships between 21 binary 
network diagnostics: degree assortativity (a), geodesic node 
betweenness (B), closeness centrality (Cc), clustering coeffi- 
cient (C), communicability (Co), geodesic edge betweenness 
(Be), eigenvector centrality (Ce), global efficiency (Eg), hi- 
erarchy (/i), local efficiency (Ei), modularity optimized using 
the Louvain (Ql) and spectral (Qs) heuristics, number of 
nodes (N), number of edges (D), mean connection distance 
(mcd), random- walk node betweenness (Brw), Rent exponent 
(p), robustness to random attack (Rr), robustness to targeted 
attack (Rt), subgraph centrality (Cs), and synchronizability 
(s). We order the diagnostics to maximize the correlation 
along the diagonal for better visualization of highly corre- 
lated groups of diagnostics. Color indicates the correlation 
between global network diagnostics over the 17 experimen- 
tal runs. (B) Reliability, as measured by the coefficient of 
variation (CV), over the 17 runs for the 21 binary network 
diagnostics reported in (A). 



APPENDIX C: COMPARISON OF CONTACT 
NETWORKS TO RANDOM GEOMETRIC 
GRAPHS 

Many of the diagnostics that we compute for the granu- 
lar networks are highly correlated with one another (see 
Fig. |9|^). In the (binary) contact networks, they form 
roughly two families, where the correlations among the 
diagnostics in a given family are high. The diagnos- 
tics that we used for the weighted networks also exhibit 
some correlations (see Fig.jsJ^), and (unsurprisingly) this 
is particularly evident for diagnostics that are known 
to be closely related mathematically. For example, the 
weighted clustering coefficient is highly correlated with 



16 



Binary Contact Network Diagnostic 


Variable 


CV 


Mean Connection Distance 


mcd 


0.0046 


Geodesic Edge Betweenness 


Be 


0.0090 


Global Efficiency 




0.0114 


Rent Exponent 


P 


0.0171 


Random- Walk Node Betweenness 


Brw 


0.0176 


Number of Nodes 


N 


0.0179 


Geodesic Node Betweenness 


B 


0.0213 


Modularity: Louvain Optimization 


Ql 


0.0071 


Modularity: Spectral Optimization 


Qs 


0.0252 


Closeness Centrality 


Cc 


0.0267 


Number of Edges 


D 


0.0360 


Synchronizability 


s 


0.0408 


Robustness, Random 


Rr 


0.0456 


Clustering Coefficient 


C 


0.0496 


Subgraph Centrality 


Cs 


0.0528 


Local Efficiency 


El 


0.0573 


Robustness, Targeted 


Rt 


0.0727 


Communicability 


Co 


0.0829 


Hierarchy 


h 


0.1080 


Eigenvector Centrality 


Ce 


0.1150 


Degree Assortativity 


a 


0.3219 



TABLE I. Reliability of network diagnostics for the binary 
contact networks. We measure reliability using coefficient of 
variance (CV). 



Force- Weighted Contact Network Diagnostic 


Variable 


CV 


Transitivity 


T 


0.0339 


Clustering Coefficient 


Cw 


0.0549 


Geodesic Node Betweenness 


Byj 


0.0282 


Geodesic Edge Betweenness 


Bew 


0.0199 


Normalized Strength 


S' 


0.0000 


Modularity: Louvain Optimization 


Qw 


0.0135 


Diversity 


V 


0.0412 


Global Efficiency 


Ew 


0.1094 



TABLE II. Reliability of network diagnostics for the force- 
weighted contact networks. We measure reliability using co- 
efficient of variance (CV). 



transitivity, and geodesic node betweenness is highly cor- 
related with geodesic edge betweenness. 

It is important to think about the possible origins 
of correlations between the various network diagnostics. 
Some of them might be specific features of the precise 
granular system under consideration, but others might 
arise because our granular packings are confined in 2D 
rather than in 3D or because of our particular prepara- 
tion protocol. Still others might be general properties of 
spatially-embedded systems (in any dimension), of gran- 
ular materials, or of networks in general. 

To examine such issues, it is desirable to compare net- 
work diagnostics computed for the networks obtained 
from experimental data with those obtained from appro- 
priate ensembles of null- model networks. It is common to 
compare the structures of networks under study to those 
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B 
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Pi 

0.1 




weighted network variables 



FIG. 9. [Color online] (A) Relationships between 8 di- 
agnostics applied to the weighted networks: geodesic node 
betweenness (Bw), clustering coefficient (Cw), diversity (V), 
geodesic edge betweenness (Bew), global efficiency (Ew), mod- 
ularity (Qw) optimized using the Louvain method, and nor- 
malized strength {S'). We ordered the diagnostics to maxi- 
mize the correlation along the diagonal for better visualization 
of highly correlated groups. Color indicates the correlation 
between global network diagnostics over the 17 experimen- 
tal runs. (B) Reliability, as measured by the coefficient of 
variation (CV), over the 17 runs for the 8 weighted graph 
diagnostics reported in (A). 



that would be expected in Erdos-Renyi random graphs 
or from some other random graph ensemble [l]. The 
networks that we study in the present paper — namely, 
contact networks in granular packings — are spatially 
embedded (in the plane) because of physical constraints. 
The development of null-models is a wide open area of 
research for spatially-embedded networks [8 , but we can 
make some progress for the binary contact networks by 
comparing the network diagnostics in those networks to 
computations of the same diagnostics using an ensemble 
of random geometric graphs (RGGs). As we will now 
discuss, we find that all diagnostics (except for the ones 
that we fix when defining the RGG ensemble to match 
their counterparts in the real networks) are significantly 
different in the real versus random networks. 

The simplest RGG |8l[32l[92] contains N nodes that are 
randomly and distributed according to some probability 
distribution throughout an ambient space, which in our 
case is R [94^ . One then places an edge between any pair 
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Diagnostic Name 






t 


P 


EG Less Than RGG 










Local Efficiency 


0.33 


0.67 


73.33 


3.1 X 10"^' 


Modularity: Louvain Optimization 


0.81 


0.92 


70.52 


7.2 X 10"^^ 


Clustering Coefficient 


0.26 


0.54 


85.28 


2.5 X 10"^° 


Random- Walk Betweenness 


0.09 


0.20 


39.73 


8.2 X 10 


Degree Assortativity 


n 1 /I 

U.14 


U.o ( 


OO.OZ 


z.iy X lu 


Subgraph Centrality 






zz. 41 


o.y X ID 


Communicability 


0.15 


0.73 


17.79 


3.0 X 10 


Rent Exponent 


0.47 


0.50 


15.41 


I.I X lU 


EG Greater Than RGG 










Global Efficiency 


0.10 


0.03 


94.41 


1.0 X 10-^^ 


Mean Connection Distance 


39.60 


32.32 


79.38 


2.5 X 10"^^ 


Robustness, Random 


8.38 X 10^ 


2.33 X 10^ 


54.89 


3.0 X 10"^^ 


Closeness Centrality 


7.84 


4.47 


50.29 


4.9 X 10"^^ 


Synchronizability 


0.0014 


0.0004 


44.71 


2.2 X 10"^° 


Robustness, Targeted 


7.08 X 10^ 


1.96 X 10^ 


36.96 


7.9 X 10"^^ 


Edge Betweenness 


13.33 


4.55 


34.71 


5.6 X 10"^'^ 


Geodesic Node Betweenness 


6.24 X 10^ 


2.16 X 10^ 


30.95 


2.0 X 10"^^ 


Hierarchy 


0.76 


0.27 


24.42 


2.9 X 10"^^ 


Eigenvector Centrality 


0.0172 


0.0080 


19.15 


4.2 X 10"^^ 


Modularity: Spectral Optimization) 


0.78 


0.75 


5.06 


1.63 X 10"^ 



TABLE in. Comparison of binary network diagnostics in the (real) experimental graphs (EGs) and the random geometric 
graphs (RGGs). We show the mean values of the EGs (column 1), the mean values of the RGGs (column 2), the t- values 
(column 3), and p- values (column 4) for a two-sample t-test between the network-diagnostic values of the two families of 
networks (EGs and RGGs). 

of nodes i and j that are separated by a distance of at 
most 2r, where one should think of the parameter r as 
the radius of a ball (using some metric) in the confining 
space. In planar Euclidean space, one considers a disk in 
and uses ordinary (Euclidean) distance. 

To compare the networks that we study to RGGs, we 
generated RGGs in which we placed nodes randomly and 
uniformly within the 2D space of the granular packing. 
For each experimental run, we created an ensemble of 100 
RGGs in which the number of nodes was identical to that 
in the experimental system. We likewise fix the number 
of edges in each RGG to be identical to that in the real 
system {D) by choosing the threshold 2r so that the num- 
ber of inter- node distances less than 2r is equal to D. We 
calculate the other 19 binary diagnostics (i.e., except for 
the number of nodes and the number of edges, as those 
have been fixed to be equal in the two sets of networks) 
and computed their mean over the 100 RGGs in each 
ensemble. By performing these computations for each 
experimental run, we created 1 estimate of each of the 
19 diagnostic values for each of the 17 runs. We report 
in Table IIIII the mean values for both the real networks 
and the networks generated from the RGG ensembles. 
We also report t-values and p-values for two-sample t- 
tests between the values in the 17 real networks and the 
17 mean values in the RGG networks. As we show in 
Table |III[ each of the 19 network diagnostics is signif- 
icantly different between the two groups. Measures of 



local connectivity (e.g., clustering coefficient) are higher 
in the RGG, whereas measures of global connectivity and 
physical distance are lower. These results illustrate that 
the networks in the RGG ensemble have more locally con- 
nectivity structures than those in the 2D granular system 
under study. 
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